Dynamics of a Josephson Array in a Resonant Cavity 
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We derive dynamical equations for a Josephson array coupled to a resonant cavity by applying 
the Heisenberg equations of motion to a model Hamiltonian described by us earlier [Phys. Rev. 
B 63, 144522 (2001); Phys. Rev. B 64, 179902 (E)]. By means of a canonical transformation, 
we also show that, in the absence of an applied current and dissipation, our model reduces to 
one described by Shnirman et al [Phys. Rev. Lett. 79, 2371 (1997)] for coupled qubits, and that 
it corresponds to a capacitive coupling between the array and the cavity mode. From extensive 
numerical solutions of the model in one dimension, we find that the array locks into a coherent, 
periodic state above a critical number of active junctions, that the current-voltage characteristics of 
the array have self-induced resonant steps (SIRS's), that when Na active junctions are synchronized 
on a SIRS, the energy emitted into the resonant cavity is quadratic in Na, and that when a fixed 
number of junctions is biased on a SIRS, the energy is linear in the input power. All these results 
are in agreement with recent experiments. By choosing the initial conditions carefully, we can drive 
the array into any of a variety of different integer SIRS's. We tentatively identify terms in the 
equations of motion which give rise to both the SIRS's and the coherence threshold. We also find 
higher-order integer SIRS's and fractional SIRS's in some simulations. We conclude that a resonant 
cavity can produce threshold behavior and SIRS's even in a one-dimensional array with appropriate 
experimental parameters, and that the experimental data, including the coherent emission, can be 
understood from classical equations of motion. 

PACS numbers: 05.45.Xt, 79.50.-fr, 05.45.-a, 74.40. +k 



I. INTRODUCTION. 

A iang-standing goal of experimentalEl^i and theoreti- 
calETEj research on Josephson junction arrays has been to 
develop sources of coherent microwave radiation. The ba- 
sic idea underlying this work is that a Josephson junction 
is a simple way of converting a d. c. current into an a. c. 
voltage. Thus, an array of A^ Josephson junctions oscil- 
lating in phase should produce a signal with N times the 
voltage amplitude, and hence A^^ times the emitted a. c. 
power, of a single junction. Arrays of overdamped junc- 
tions have seemed most promising for coherent emission, 
since junctions of this type have, at any given applied 
current, only a single voltage state, and thus have none 
of the multistability and chaotic behavior which could 
inhibit coherent emission. However in practice, it has 
proved very difficult to achieve an efficient a. c. to d. c. 
conversion in such systems; thus-jfac, the highest conver- 
sion efficiency is only about l%Ej^E3. The low efficiency 
may result from the high degree of neutral stability which 
has been shown to exist in such overdamped arrays in the 
absence of an applied magnetic fieldla. 

Recently, a remarkably high conversion efficiency has 
been achieved in an underdamped Josephson array, by 
coupling the junctions in that array electraaiagnetically 
to a mode in a resonant microwave cavityllll. The high 
emission is a manifestation of the so-called self-induced 
resonant steps (SIRS's) that appear on the current- 
voltage (IV) characteristics of these arrays. It is thought 
such arrays emit strongly because every junction is cou- 
pled to the same electromagnetic mode, and hence, ef- 



fectively to every other junction. This same coupling 
is presumed to lead to the observed threshold effect, in 
which the strong emission occurs only above a certain 
array length. 

A number of models have been proposed which pro- 
duce some aspects of this behavior. In the original model, 
which inspired the measurements, an analogy was drawn 
between junctions in a voltage-biased series Josephson ar- 
ray and a collection of two-level atoms where eepulation 
inversion and laser emission could be achievedtS. Several 
authors have introduced various kinds of impedance loads 
across groups of junctions or a one-dimensional array, in 
order to achieve a global coupling and hence, to,- 
gate coherence among the junctions in the arrayll3ll^ 
None of these models have yet produced both the self- 
induced resonant steps and the threshold junction num- 
ber seen in the experiments. The coupling of prop- 
agating modes in Josephson ladders and other struc- 
tures to electromagnetic radiation has also been stud- 
ied theoreticallyn3Ej. A simple Hamiltonian to treat 
the equilibrium properties of a one-dimensional, voltage- 
driven array in the weak-coupling regime has recently 
been proposed and studied within a mean-field approxi- 
mation which should be, very accurate in the limit of large 
numbers of junctionjcJ. Within this approach, it was 
found that the array developed coherence only above a 
threshold number of juni±ions, in agreement with exper- 
iment. In a recent papeiH, the present authors proposed 
a similar model to treat array dynamics for any strength 
of coupling; they also briefly described a few numerical 
results obtained from the model, including both a thrcsh- 
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FIG. 1. Sketch of the array geometry considered in our 
model. There are junctions (crosses), labeled by their 
gauge- invariant phase differences 7^, and A'' + 1 supercon- 
ducting islands. A current / is injected into one end of the 
array and extracted from the other. The array is placed in 
an electromagnetic cavity which supports a single resonant 
photon mode of frequency Q. 

old for coherence and self-induced resonant steps. 

In the present paper, we give a more complete deriva- 
tion of the model equations of motion of Ref. [25] for a 
one-dimensional array of underdamped Josephson junc- 
tions coupled to a single-mode electromagnetic cavity. 
Starting from a suitable Hamiltonian, we obtain the 
Heisenberg equations of motions for the phase differ- 
ences and the photon creation and annihilation opera- 
tors. We account for dissipation in the junctions by the 
standard procedure of coupling each junction to a reser- 
voir of phonon variables with a density of states con- 
structed so as to produce ohmic damping. In the limit of 
large numbers of photons, the equations can be treated 
classically and solved numerically. We also correct the 
treatment of Ref. [25] of the junction damping and the 
coupling of the array to an external current. Finally, we 
carry out a canonical transformation of our Hamiltonian 
to show that the interaction between the array and the 
cavity mode has the form of a capacitive coupling. 

We also present much more extensive numerical results 
than those of Ref. [25], based on solutions to the model 
equations. Our numerical results show all the princi- 
pal features of the measurements, including SIRS's, a 
coherence threshold, and a quadratic dependence of the 
photon energy in the cavity upon number of active junc- 
tions. Plots of the IV characteristics and other calculated 
features closely resemble the corresponding experimental 
plots. This agreement is especially noticeable since the 
calculation is one-dimensional, while most experiments 
have been conducted on two-dimensional arrays. In ad- 
dition, we find that by tuning the initial conditions, we 
can cause the array to lock into a variety of different in- 
teger SIRS's, again in agreement with experiment. We 
conclude that these equations do indeed describe the ex- 
periment, and that a one-dimensional array is sufficient 
to achieve this type of coherent behavior. 

The remainder of this paper is organized as follows. In 
Section II, we derive the Heisenberg equations of motion 



for the phase and photon variables, starting from a model 
Hamiltonian. We also apply a canonical transformation 
which shows that the Hamiltonian involves a kind of dis- 
tributed capacitive coupling between the Josephson array 
and the cavity mode. In Section III, we give a detailed 
description of our numerical results. Section IV presents 
a comparison between our results and experiment, gives a 
qualitative discussion of the numerical results, and makes 
some concluding remarks about the model. 

II. DERIVATION OF THE EQUATIONS OF 
MOTION 

A. Hamiltonian 

We consider a one-dimensional array of N Josephson 
junctions placed in a resonant cavity, which we assume 
supports only a single photon mode of frequency fl (the 
geometry is sketched in Fig. |lj). The array is to be driven 
by an applied current /. We write the Hamiltonian in the 
form 

H = H photon + Hj + He + Hcurr + Hdiss- (1) 

Here, Hphoton is the energy of the cavity mode, which we 
express as 

Hphoton = Tt-VI (^^a + , (2) 

with and a as the usual photon creation and annihi- 
lation operators. Hj is the Josephson Hamiltonian, and 
is assumed to take the form 

N 

Hj = -Y,Ejj COS jj, (3) 
i=i 

where Ejj is the Josephson energy of the j*'* junction, 
and 7j is the gauge-invariant phase difference across the 
j*'* junction (defined more precisely below). Ejj is re- 
lated to Icj, the critical current of the j*'* junction, by 
Ejj ~ hicj/q, where q = 2|e| is the Cooper pair charge. 
He is the capacitive energy of the N junctions, which we 
approximate as 

N 

Hc^Y.^o,nl (4) 

Here, Ecj = q^ j {^Cj) is the capacitive energy of the j*'' 
junction, Cj is the capacitance of that junction, and Uj 
is the difference between the number of Cooper pairs on 
the j*'' and (j -I- 1)'^ superconducting island. 

The gauge- invariant phase difference, 7-,-, is the term 
which couples the Josephson junctions to the cavity. It 
may be written as 
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7, - 0, - [(27r)/<i>o] / A • ds ^ 



(5) 



prevent the original potential from being sliifted by the 
coupling to the j*'* phase degree of freedorrH. The spec- 
tral density of the harmonic oscillators in the j*^ junction, 



where (jij is the phase difference across the j*'' junction denoted J]{uj), is defined by 



in a particular gauge, A is the vector potential in that 
same gauge, $o = hc/q is the flux quantum, and the line 
integral is taken across the junction. We assume that A 
arises from the electromagnetic field of the normal mode 
of the cavity. , \^ [Gaussian units, this vector potential 
takes the forr 



A(x, t) = ^{hc^)l{n) {a{t) + a\t)) E(x), 



(6) 



where E(x) is a vector proportional to the local electric 
field of the mode, normalized such that Jy (i^x|E(x)p = 
f , where V is the cavity volume. Given this representa- 
tion for A, the phase factor Aj can be written 



where gj takes the form 



9j 



m (27r) 
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E(x) • ds 



(7) 



(8) 



nthe 



Clearly, gj is an effective coupling constant describiit 
interaction between the j*'' junction and the cavityE3. 

The terms discussed so far need to be supplemented 
by the effects of a driving current and of damping within 
the junctions. A driving current is easily included within 
the Hamiltonian formalism via a "washboard potential," 
Hcwrr, of the form 



hi 



N 



(9) 



i=l 



with / as the driving current. 

Tiie inclusion of dissipation can be done in a standard 
wayE3~E3 by coupling each gauge-invariant phase differ- 
ence, 7j, to a separate collection of harmonic oscillators 
with a suitable spectral density. Thus, we write the dis- 
sipative term in the Hamiltonian as 



N 



(10) 



where 



n2 . 
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U 



'a,j 



2m 



-(7, 



(11) 



The variables Ua.j and Pa.j, describing the a oscillator 
in the j*'* junction, are canonically conjugate, and fUaj 
and LOaj are the mass and frequency of that oscillator. 



The last term in Eq. (11) must be added in order to 



maUJa 



S{uj — LOa)- 



(12) 



If Jj{i^) is linear in \uj\^ it can. he. shown that the dissipa- 
tion in the junction is ohmicEjLJ. We write such a linear 
spectral density as 



2ti 



\uj\ e(wc - w), 



(13) 



where ujc is a high-frequency cutoff (at which the as- 
sumption of ohmic dissipation begins to break down), 
@{u>c — Lo) is the usual step function, and aj is a dimen- 
sionless constant, which we write as aj — Ro/Rj, where 
Ro = ft-/(4e^) and Rj is a constant with dimensions of 
resistance (actually, the effective shunt resistance of the 
junction, as discussed below). 



B. Equations of Motion 

It is now convenient to introduce the operators and 
ai by 



a = ajf + iaj] 



aji - laj. 



(14) 



(15) 



The free photon part of the Hamiltonian can be expressed 
in terms of aj^ and aj as follows: 



Hphoton — flfl{(lji 



(16) 



where we have used the additional commutation rela- 
tions [ail, o/] — */2i which follows from the usual relation 
[a, a^] = 1. The gauge-invariant phase difference, jj, is 
related to 4>j by 



Ij = '/'i - '^\/TjaR. 



(17) 



The time-dependence of the various operators appear- 
ing in the Hamiltonian (|l|) is now obtained from the 
Heisenberg equations of motion. For a general operator 
O, these take the form 



(18) 



These equations of motion can be evaluated for the var- 
ious operators entering H, using the commutation rela- 
tion [A, F{B)] = [A, B]F'{B), where F is any function of 
an operator B, and F'{B) is the derivative of that func- 
tion. One also needs the commutation relations for the 
various operators in the Hamiltonian (^. Besides the 
relations already given, these are as follows: 
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[Paj-,U/3,fc] = -i^ 5a,P Sj^k- 



(19) 
(20) 



Note that 7^, unlike no longer commutes with a/; 
instead, it satisfies 



[lj,aR\ = 0; 



/9j- 



(21) 
(22) 



Using all these relations, we find, after a little algebra, 
the following equations of motion for the operators jj, 
rij, an, and a/: 



7, =2 



Ecj_ 
h 



a/. 



sin(7j) + - 
n q 



a \ 



(/a 



a/ = -f^ fli? + 51 V57 sin(7j ) - - ^ 



/9j 



2 



(23) 

(24) 
(25) 

(26) 



These are equations of motion for the operators or, ai, 
rij, and (or 7^). Note that they do not depend on 
the particular choice of gauge, but only on the form of 
the Hamiltonian and the commutation relations for the 
various operators. We will study these general equations 
within the limit of large number of photons in the cav- 
ity and large number of charges in the junctions, and 
in this "classical" limit, we will regard the operators as 
c-numbersE2l. 

We also have the equations of motion for the harmonic 
oscillator variables. Since we have no explicit interest 
in these variables for themselves, we instead eliminate 
them in order to incorporate the dissipative term into 
the equations of motion. Such a replacement is possible 
provided that the spectral density of each junctien-^ lin- 
ear in frequency, as noted above. In that caseE3"l23, the 
oscillator variables can be integrated out. The effect of 
carrying out this procedure is that one should make the 
replacement 



wherever this sum appears in the equations of motion. 

Making the replacement (|7|) in Eqs. (H) and (|26|), 
we obtain the equations of motion for rij and a/ with 
damping: 



h q 2uJcjQjj 



7j 



(28) 



ai - 



3 



3 



(29) 



Here, we have introduced the parameters ojcj = E^j/h, 
which is a frequency associated with the capacitive en- 
ergy of the j*'' junction; Cup = ^ SjLi ^pj^ the average of 
the Josephson plasma frequencies ujpj — ^2EcjEjj/fi] 
and Qjj, the dimensionless junction quality factor (or 
damping parameter) for the j*'* junction, which is related 
to the capacitance Cj and the shunt resistance Rj by 



LOpRjCj. 



(30) 



Eqs. dH), (HI), (^), and (^ ) can be combined, with a 
little algebra, into two coupled second-order differential 
equations: 



2uj. 



2ujcjQjj 



and 




ujj sm7j 

L _ • 



aR 



(31) 



"'«^ = -i E^^- (32) 



UJCj 



where we have defined ujjj — Ejj/h. Note that in the 
absence of coupling between the junctions and the cavity, 
7i ~ '^ii -^1- (0) reduces to the usual equation of m£t 
tion for a resistively and capacitively shunted junctionEj 
driven by a current /, and Eq. ( ^2| ) reduces to that of 
a simple harmonic oscillator which represents the cavity 
mode. Note that we have not included any damping due 
to the cavity walls. While such damping is undoubtedly 
present, we find that good agreement with experiment 
can be obtained without including it. 



C. Canonical Transformation 

The physics behind the coupling between the Joseph- 
son junctions and the resonant cavity, and hence the 
physics of Eqs. (^ ) - (^, can be made clearer by a 
canonical transformation. For simplicity, we describe 
this transformation including only the terms Uphot, Hj, 
and He from the Hamiltonian (|l|), and omitting Hcurr 
and Hdiss ■ The same transformation has previously been 
used for a single junction coupled to a resonant cavity by 
Buisson and HekkingEJ; and, for two voltage-driven junc- 
tions coupled to a resonant cavity, by Shnirman et aO. 

We begin by writing 
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H' = Hphot+Hj + Hc = -hn{pl + q^.) 



N 

E 



Ecjfi^j - Ejj cos{(f>j 



(33) 



where we have defined pr = o/ and qr = an- With 
this choice, Pr and qr satisfy the commutation relation 
[PriQr] = —i- Next, we make the canonical transforma- 
tion 



(34) 
(35) 




Qr = Qr- 



(36) 
(37) 



The only nonvanishing commutators of the primed vari- 
ables are easily shown to be \n^, (p'^] — [pj,, q^] = —i. Re- 
expressing the Hamiltonian ( p3|) in terms of the primed 
variables, we obtain 



N 

E 



Ec]{n'jf - Ejj cos0^- 



(38) 



Thus, H' is the sum of four terms: the sum 
{hQ/2)[{pr)'^ + {q^)^] describes the cavity resonator; the 
last sum describes the N independent junctions; and the 
remaining terms represent the interaction between the 
junctions and the resonator, and an indirect interaction 
between the junction variables n'^ mediated by the cavity. 

To interpret this interaction, we note that the junction- 
cavity system has two places in which charge can be 
stored: the variables n'^ of the junctions and the vari- 
ables p'j. of the cavity. The cavity behaves much like an 
LC circuit, with capacitive energy (?iil/2)(pj,)^ and in- 
ductive energy (?ir2/2)(g^)^. The dominant interaction is 
a capacitive coupling between the charge variables nLs)f 
the junction and the charge variable pj, of the cavityE3. 

In further support of this interpretation, we now show 
that H' , in the-iorm (38), is equivalent to that given by 
Shnirman et aSB in the case of zero applied voltage. Fig. 
1 of Ref. [35] depicts each junction contained between the 
plates of a capacitor with capacitance C. The junction 
itself has capacitance Cj. The system of junction and 
capacitor is then shunted by a parallel inductance, L, 
which acts to couple all the junctions togetheiO. The 
equivalence is established by noting the correspondence 
between the variables used in the present paper and the 
variables 0, q, rij and 9j of Ref. [35]. The correspondence 
is as follows (assuming that the coupling constants gj, 
Ecj, and Ejj are independent of j): 



q'r ^ [L/(2C0]-^/V; 

p; ^ [L/{2Ct)Y'\, 



Hi 



where Ct = CCj/{C + 2Cj). In order to complete 
the correspondence, we also give the correspondence be- 
tween the remaining parameters and the quantities L, C, 
Ct,and Cj: 



2g- 
n = 

Ec^ 



L 



1/4 



2CtJ Cj' 

1 



V2LCt' 
1 

C + 2C,j 



If we make the replacements and identifications given 
above, then our Hamiltonian, for the case of two junc- 
tions, is identical to that considered in [35] in the absence 
of an applied voltage. The main differences between the 
two models are the boundary conditions: constant cur- 
rent bias in our model, and fixed voltage bias in that of 
Ref. [35]. 



D. Dimensionless Form 

In order to write these equations of motion in a simpler 
form, we now introduce the dimensionless time t — Qpt. 
Also, although we allow disorder in the parameters of 
the different Josephson junctions, we assume that the 
coupling constants between the junctions and the cavity 
are all the same, i. e., that gj — g for all j. In addition, 
for simplicity, we assume that the products RjCj and 
Icj/C'j, and hence ujpj, are independent of j. Then Qjj 
is also j-independent, and we may introduce the dimen- 
sionless coupling constant 



5 = 5 



3l 

TlUJn 



(39) 



which is also independent of j. Similarly, we introduce a 
scaled number variable fij by 



2Ec^ 

TlUJr, 



a dimensionless frequency f2 by 
photon variables and a/ by 

= V5 an; 
ai ^ y/g aj. 



(40) 

il/ujp, and scaled 

(41) 
(42) 



Finally, we assume that the critical currents Icj are ran- 
dom and uniformly distributed between /c(l ~ A) and 
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FIG. 2. Left-hand scale and solid line: current-voltage (IV) 
characteristics of a one dimensional array of A*' = 40 junctions, 
with disorder parameter A = 0.05 and coupling constant 
g — 0.001. The resonant frequency of the cavity is f2 = 2.2, 
and the damping parameter of the junctions is Qj = \/20. 
Right-hand scale and stars : scaled total energy E = gE/{fiQ.) 
carried by the resonant mode of the cavity, plotted as a func- 
tion of decreasing current I /Ic- The vertical dashed lines 
are guides to the eye. The upper dashed horizontal line indi- 
cate the expected position of the integer self-induced resonant 
steps (SIRS's) for the particular resonant frequency of the 
cavity (all junctions in the n = 1 SIRS). For the lower dashed 
horizontal line, 23 junctions are on the n = 1/2 SIRS, and 17 
junctions are in the = state. Branches correspond- 

ing to increasing and decreasing current are shown by arrows. 
Double-headed arrows on this and subsequent Figures denote 
that the curve can be obtained by sweeping the current in 
either direction. 

/c(l + A), where A is a measure of the degree of disorder. 
After some algebra, we eventually obtain the following 
equations of motion: 



I 



' Ic{l- 



— ^ - sin(7j) + 2 — a/; 



Qj 



Qj 



(43) 
(44) 
(45) 



aj = —ft aji — 2 ft g — ^ 
Qj 



Ed 



N g 



+ A,) sin(7,) + + A,) n,. (46) 



In these equations, the dot refers to differentiation with 
respect to r, and the j*'* critical current is /c(l + Aj). 

These equations can be combined into two more com- 
pact equations, with the result 



2 an; 



(47) 



+ {n'f ~aR = -~g-^ ^(1 + A,) 7,, (48) 

where we have defined (f7')^ fiV[l+2 5 

Eqs. (H) and (||) are the analogs of Eqs. ^ and (|^), 

expressed in terms of dimensionless reduced variables. 



III. NUMERICAL RESULTS 



We have solved Eqs. (47) and (W8[) for the variables 



hj, 7j, qr and a/ numerically using the same approach 
as in Ref. [25], namely, by implementing_.the rapid and 
accurate adaptive Bulrisch-Stoer methocO. We initialize 
the simulations with all the phases randomized between 
[0, 27r], and usually an ~ aj = hj = 0. We then let the 
system equilibrate for a time interval Ar = 10"*, after 
which we evaluate averages over a time interval Ar = 
2 • 10'^, using 2^^ evenly spaced sampling points. 



A. Typical IV Characteristics, Power Spectrum, and 
Coherence Transition 

In Fig. ^ we show a representative current-voltage (IV) 
characteristic calculated for an array of TV = 40 junctions 
with A = 0.05 and g = 0.001. The time-averaged voltage 
{V)r [left-hand scale] is obtained from 



N 



(49) 



where {■■■)t denotes a time average and Vj is obtained 
from the Josephson relation. 



ii. 

RIc 



h d'jj 1 



(50) 



A striking feature of this plot is the self-induced reso- 
nant steps (SIRS), at which {V)t remains approximately 
constant over a range of applied current. For this par- 
ticular choice of parameters and initial conditions, we 
see these steps at {V) r / [N RIc) = nVl/Qj. These steps 
corresponds to voltages at which the condition 



2e{V.j)t = nUn 



(51) 



is satisfied for the individual junctions, with n — 1 (up- 
per horizontal dashed line) and n = 1/2 (lower horizontal 
dashed line). Thus, the lower step is at 23/80 the volt- 
age of the upper step. For the latter case, the driving 
current is smaller than the retrapping currents of 17 of 
the junctions; thus, only 23 out of the 40 junctions are 
oscillating on this step. (The retrapping current is the 
minimum current for which an underdamped junction is 
bistable.) The steps occur at exactly the voltages where 
the first integer and half-integer steps would appear in 
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FIG. 3. Power spectrum, P{uj) (Eq. of the a. c. voltage 
across the array, plotted versus frequency, fl, at two driving 
currents: (a) and (b) I/Ic = 0.58, corresponding to the first 
integer SIRS, and (c) and (d) Ijl^ = 0.65^ slightly off a SIRS. 
Other parameters are the same as in Fig. pi Panels (b) and (d) 
are the same as (a) and (c) except that the effective coupling 
to the resonant cavity, g = 0. In each panel, the left verti- 
cal dashed line shows the resonant frequency of the cavity, 
and the right vertical, dashed line shows the average resonant 
frequency of the junctions for the case of no coupling to the 
cavity. 

these junctions, if the junctions were driven by an a. c. 
current of frequency J7. Thus, the radiated energy in the 
cavity seems to behave hke an a. c. drive which acts back 
to induce these steps in the junctions of the array. Sim- 
ilar steps were seen experimentally in a two-dimensional 
array of underdapped Josephson junctions coupled to a 
resonant cavityll3, and in more recent experiments in ID 
arraysO. 

Fig. H also shows the time-averaged scaled total en- 
ergy, E, contained in the cavity, [right-hand scale of the 
Figure]. E is defined as 



E = {a-n + aj)r = g (a^ + aj)r 



_9_ 



E, 



(52) 



where E = {Hphotonjr is the cavity energy; it is plotted 
as a function of 1 1 Ic for the same array. As is evident, E 
increases dramatically when the array is on a SIRS, and 
is very small otherwise. This sharp increase signals the 
onset of coherence within the array, and can be qualita- 
tively understood from the equations of motion. Specifi- 
cally, when the array sits on one of the integer SIRS, all 
the junctions are oscillating in phase. Hence, the term 
driving or [the right-hand side of Eq. (^)], and thus aji 
itself, are both proportional to the number Na of active 
junctions. 

Before proceeding further, we briefly review the con- 
cept of active junction number Na, as discussed in Refs. 
[12] and [25]. This concept has meaning only for under- 
damped junctions. Such a junction is bistable and hys- 



teretic in certain ranges of current - that is, it can have 
either zero or a finite time-averaged voltage across it, de- 
pending on the initial conditions. In the present case, Na 
denotes the number of junctions (out of N total) which 
have a finite time-averaged voltage drop. It is possible 
to tune Na by suitabl y (j ho osing the initial conditions, 7i 
and 7i, in simulationsPtri. 

Fig. ^ shows the calculated voltage power spectrum of 
the a. c. component of the total voltage across the array 



P(c.) 



2 lim 



(53) 



for two values of the driving current: I/Ic = 0.58 [Fig. 
^(a) and (b)] and I/I^ = 0.65 [Fig. |(c) and (d)]; aU other 
parameters are the same as in Fig. ^. In (a), all the junc- 
tions are on the first SIRS, while in (c), the array is tuned 
off this step. In Fig. ||(b) and ^(d), we show the same 
case as in Fig. ||(a) and (c) respectively, except that the 
coupling constant, g, is artificially set equal to zero. Note 
that in ||(a), the power spectrum has peaks at the scaled 
cavity frequency, fi, and its harmonics. This is evidence 
that the junctions are all oscillating at frequency fl. In 
case (b), the junctions are still coupled by the indirect 
interaction via the cavity, but the power spectrum shows 
that the array is not synchronized in this case; instead, 
the individual junctions oscillate approximately at their 
individual resonant frequencies and their harmonics and 
subharmonics. Hence, the power spectrum has a spread 
of frequencies, all of which differ from that of the cavity. 
In cases (b) and (d), the junctions are, of course, inde- 
pendent of one another, and the power spectrum is that 
of a disordered one-dimensional Josephson array with no 
coupling between the junctions. 

We have also calculated the response of a disordered 
array (A = 0.05) of fixed length {N = 40 junctions), and 
a driving current I/Ic = 0.58, when the number of active 
junctions, Na is varied. This current not only lies well 
within the bistable region, but also leads to a voltage 
on the first integer SIRS. In Fig. §(a), we plot the time- 
averaged scaled energy of the cavity, E{Na) [Eq. (^)], as 
a function of Na- For Na < 17, the active junctions are 
unsynchronized, and E is correspondingly small and only 
weakly dependent on Na- There is a sudden jump in E 
at a critical number of active junctions Nc = 17. Above 
this value E increases as a quadratic function of Na , and 
we have fitted E{Na) to the form E = co + ciNa + C2Na- 
The constants which give the best fit are cq = —0.00163; 
ci = 0.00125; and C2 = 6.868 • 10"^. This curve is shown 
as a full line in Fig. ||(a); the fit is clearly excellent. As 
a contrast, we also show the best linear fit to the same 
data set (dashed line); the fit is plainly less good. The 
magnitude of the jump in E at Na = Nc is nearly a factor 
of ^ 10'^, as shown in the inset to figure. 

To measure the degree of synchronization among the 
Josephson junctions, we, have also calculated the Ku- 
ramoto order parametertB, {r)r, for the same parameters, 
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FIG. 4. (a) Asterisks: Scaled photon energy E — gE/{hQ) in the resonant cavity when the array is current driven on a SIRS, 
plotted versus number of active junctions, Na- The array parameters are A*' = 40, Q — 2.2, Qj — y/20, A = 0.05, g — 0.001, 
and I/h = 0.58 [cf. Fig. | (a)]. FuU curve shows the best fit of E to the function C2N^ + CiNa + Co for Na > 17, the threshold 
for synchronization. The fitting parameters are co = -0.00163, ci = 0.00125, and C2 = 6.868 • 10"''. We contrast this fit to 
the best linear fit (dashed line). Inset: E(Na) near Nc — 17, showing jump near synchronization threshold, (b) Open circles: 
Kuramoto order parameter, (r),- [Eq. ([54[)], for the same array. Dots connecting circles are guides to the eye. The sharp 
increase in {r)^ and the quadratic increase in E both begin at Nc — 17. (c) Measured a. c. power as a function of the input 
d. c. power, as obtained in Ref. [17] for a 3 x 36 array. The d. c. power is proportional to the number of active rows in their 
array, while the a. c. power is proportional to the energy E in the cavity. 



as a function of number of active junctions, Na- (r) 
defined by 



^exp(i7j)|). 



(54) 



The results are shown in Fig. ^(b). Note that (r)^ = 1 
represents perfect synchronization among the active junc- 
tions, while {r)^ — would correspond to no correla- 
tions between the different phase differences, (j>i. Just 
as for E{Na), there is an abrupt increase in (r)r at 
Na — Nc, indicative of a dynamical transition from an 
unsynchronized to a synchronized state (with all active 
junctions locked to the same frequency and having a com- 
mon phase) , as Na is increased keeping all other parame- 



ters fixed. As with similar transitions in other modelscil, 
this transition is not inhibited by the finite disorder in the 
Ic's. Instead, (r)T- approaches unity, representing perfect 
synchronization, (r)^ remains finite even for Na < Nc, 
because even in this regime there is still some residual 
correlation among the phases in different active junc- 
tions. This transition is the dynamic analog of that an- 
alyzed by an equilibrium mean- field theory in Ref. [24] . 

Finally, in Fig. ^c), we show an experimental plot of 
the detected a. c. power as a function p£ the input d. 
c. power, as measured by Barbara et atll for a 3 x 36 
array. These quantities are, of course, not equivalent to 
the calculated results which are plotted in Fig. I (a). The 
input d. c. power is equal to the power dissipated in the 
active junctions; so it is proportional to Na- The detec- 
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FIG. 5. Same as Fig. |(a), except that TV = 80. In this 
case, the synchronization threshold is Nc — 20, and the 
quadratic fit to the energy above synchronization has differ- 
ent fitting parameters: Co = —0.01576, ci = 0.001149, and 
C2 = 1.441 • IQ-^ 

ted a. c. power is that measured by a pickup junction in 
the cavity, and thus should be proportional to E{Na) in 
our notation. Despite the differences, our calculated plot 
(for a one-dimensional array) appears strikingly similar 
to their measured plot, especially as regards the disconti- 
nuity at the threshold and the quadratic dependence on 
Na for Na above the threshold. 

In Fig. 1^, we show the synchronization transition for 
an array of TV = 80 junctions, keeping the other param- 
eters the same as in Fig. ^a). In this case, the criti- 
cal threshold is Nc = 20, somewhat larger than for the 
A'' = 40 junction array. The inset shows that the cavity 
energy still has a discontinuity by a factor of ~ 10"^. How- 
ever, the quadratic function which best fits E{Na) for 
Na > Nc is now described by the different fitting parame- 
ters: Co = -0.01576, ci = 0.001149, and C2 = 1.441-10-5. 
Thus, the total length of the array alters the details but 
not the qualitative features of E{Na)- 

These calculations were carried out for an array tuned 
to the first SIRS. If, instead, we carry out the same cal- 
culation when the array is tuned to the bistable region 
but not tuned to a SIRS, we find that E does not increase 
quadratically with Na- Instead, E{Na) shows no thresh- 
old behavior, and, indeed, varies little with Na- A plot 
of E{Na) in this case is shown in Fig. ||. The parameters 
are the same as for the calculation in Fig. |^(a), except 
that the driving current in this case is I /Ic = 0.65, which 
is not on a SIRS [cf. Figs. | and | ]. 



B. Effects of Varying the Number of Active 
Junctions 

In Fig. 0(a) , we show a series of IV characteristics for 
a 10-junction array {N — 10), calculated by varying the 
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FIG. 6. Total scaled cavity energy i5 as a function of the 
number Na of active junctions, for the same array parame- 
ters as in Fig. ^a) except that the current is tuned off any 
self-induced resonant step: I/Ic = 0.65 [cf. Fig. ^(c)]. In this 
case, E does not increase quadratically with Na above a criti- 
cal threshold; instead, it shows no threshold behavior, is only 
weakly dependent on Na, and is much smaller than in Fig. 
ga). 

number Na of active junctions from 1 to 10. Each solid 
vertical line segment corresponds to the IV characteristic 
for a different Na, and represents Na junctions sitting on 
the first integer SIRS. The width of each segment rep- 
resent the current height for that step, as found in our 
calculation. The dashed vertical lines show the expected 
voltages for the integer SIRS's, and are good matches for 
the calculated voltages for the various Na's. The long 
straight diagonal line segment, which is common to all 
the different Na's, represents the ohmic part of the IV 
characteristic with all junctions active. The nearly hori- 
zontal dashed line in the upper right hand corner of the 
Figure shows the IV characteristic for increasing voltage 
with Na = 10. The very short vertical segments within 
this dashed line correspond to several junctions which 
have been excited to higher steps, specifically the n = 2 
(second integer step) while the remaining junctions are on 
the n — 1 step. The horizontal dashed line on the lower 
left represents the low-voltage end of the Na = 10 IV 
characteristic (on decreasing current). The short vertical 
segment within this dashed line corresponds to fractional 
SIRS's - specifically, three of the junctions have slipped 
from the n = 1 to the n — 1/2 step, while the rest are 
in the {Vj)T = state (the driving current is smaller 
than their individual retrapping currents). Thus, we see 
both the higher integer and the fractional SIRS's in these 
one-dimensional arrays. 

In Fig. 0(a), although we show the full hysteresis loop 
only for Na = 10, the IV curves for other values of Na 
are also hysteretic. In all cases for which Nc < Na < 10, 
the number of active junctions increases when the SIRS 
becomes unstable, and individual junctions jump into the 
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FIG. 7. (a) Current I/Ic versus time-averaged voltage {V}t/{NRIc) for an array containing A'" = 10 junctions, and with 
damping parameter Qj = y/20, disorder A — 0.05, cavity coupling g — 0.003, and cavity resonance frequency of = 1.8. By 
properly choosing the initial conditions, one can select the number A^^ of active junctions to be any integer between and 10. 
Each vertical line segment corresponds to a portion of the IV characteristic for a particular choice of Na, as obtained with 
increasing current (although the same result would be obtained with decreasing current). The ohmic (straight diagonal line) 
segment is found for Na ~ 10 with decreasing current. The dashed vertical lines indicate the voltages of the expected integer 
SIRS's. The dashed, nearly horizontal line corresponds to increasing the voltage on the A^a = 10 IV characteristic; the dots and 
very short vertical line segments within this dashed line corresponds to currents at which several of the active junctions jump to 
the n — 2 SIRS. The short, nearly horizontal dashed line in the lower left-hand corner occurs on the decreasing current branch 
with 10 active junctions. The very short vertical line segments within this dashed region correspond to several active junctions 
synchronizinignon the n = 1/2 SIRS, while the remainder are in the state of {Vj)r = 0. (b). Measured IV characteristics for a 
3 X 36 arraytZI. The open circles represent self-induced-^esonant steps corresponding to different numbers of active rows. Full 
squares are believed to be examples of resistance stepsH. 



n — 2 SIRS state; ohmic behavior is not attained until 
I/Ic > 1. For Na < Nc, the array behaves somewhat 
differently: when the SIRS becomes unstable, Na is un- 
changed, and the IV curve immediately becomes ohmic. 
When I /Ic ~ 1 in this regime, the remaining junctions 
become active and the IV characteristic also becomes 
ohmic. For this particular array, Nc — 4. 

As a comparison, we also show, in Fig. |^(b), the IV 
characteristics as measured for a 3 x 36 underdamped 
array, by Barbara et aO. The open circles correspond 
to the steps observed for different numbers of active rows 
(from 1 to 10 in this instance), which are produced when 
an in-plane magnetic field reduces the critical current of 
the individual junctions. The more widely spaced daclc 
rows are believed to be examples of resistance stepsEj. 
The steps (open circles) very much resemble those of Fig. 
^(a), even including the low-current falloff (though the 
shapes of the curves are slightly different). 

We have also calculated E, the energy in the cavity, as 
a function of injected d. c. power, Pdc, when the array is 
biased on a SIRS, for several choices of array parameters. 
A typical example of our results is shown in Fig. ||(a), 
where E is plotted versus Pdc = {I / Ic)[{V)r / [N RIc)] for 
an array of ten junctions, using the same parameters as 
in Fig. ^(a) and varying the values of Na- Each curve 
corresponds to a different number Na of active junctions, 
and, for each Na, we sweep current across the n = 1 SIRS 



(leftmost curve corresponds to Na — 1, and rightmost 
to Na = 10). The curves end when the SIRS's become 
unstable. Each curve is quadratic at low Pdc and ap- 
proximately linear at higher Pdc- For commrison we also 
show the corresponding experimental plot^ for a 4 x 36 
array for Na = 16, 21, and 23 active rows [Fig. ||(b)]. In 
all cases the experimental array is above Nc, the coher- 
ence threshold. The similarity between the experimental 
and calculated curves is strikingly apparent. 



C. Effects of Changing Model Parameters 

Finally, we have studied how our numerical results de- 
pend on the parameters of our model. There are several 
parameters of interest: the number of junctions N, the 
disorder parameter A, the damping parameter Qj, the 
coupling constant g, and the normalized cavity mode fre- 
quency rj. Clearly, a thorough numerical investigation of 
all these parameters is out of the question. We have 
therefore varied only two parameters in the present pa- 
per: Qj and g. 

Fig. H shows the total time-average voltage {V)t across 
the array, and the total time-averaged energy E in the ar- 
ray as a function of driving current I/Ic, for Qj — \/100, 
\/lo, V2, VO^, and Vol35, all for g = 4 x IQ-^ N = 10, 
and A — 0.05. In each case, the resonant frequency of 
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FIG. 8. (a) Calculated total energy E within the cavity, plotted versus d. c. power, Pdc, for Na activeiunctions synchronized 
on the n = 1 SIRS, for an array of 10 junctions (A'' — 10), using the same parameters as in Fig. tA Each curve segment 
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enlargement of the calculated curve for Na ~ 6 (filled circles), (b) Experimental results for a 4 x 36 array as reported in Ref. 
[43]. From left to right, these results correspond to Na = 16, Na = 21, and A'^a = 23 active rows (all in the coherently radiating 
state with Na > Nc). 



of the cavity fl is chosen such that the scaled voltage 
^/Qj — 0.9. This choice insures that the voltage lies 
within the bistable region of the IV characteristic for the 
underdamped junctions. The arrows in the upper panel 
indicate the direction in which the current is swept. We 
show only the energy in the cavity for the decreasing 
current branch. 

Several features of these curves are apparent. First, the 
SIRS's are wider on the increasing than the decreasing 
branches. For the most underdamped case (a), there are 
no visible SIRS's on decreasing the current. Secondly, the 
cavity energy shows clear signs of a resonant interaction 
between the array and the cavity in cases (a)-(c). Finally, 
there are strong indications of an integer SIRS even for 
the overdamped case (d), where there is no bistable re- 
gion in the uncoupled IV characteristics. [We find an 
even clearer integer SIRS in case (d) if we increase g by 
a factor of 10. In this SIRS also develops in (e) 

(not shown in the Figure)]. 

In Fig. ^ (a) - (d), we plot {V)t and cavity energy E 
versus I/Ic for several values of the coupling constant g, 
all for Qj = \/20, TV = 10, A = 0.05, and n/Qj = 0.9. 
Once again, the arrows in the upper panels denote direc- 
tion of current sweep. As discussed in the next section, 
we believe that experiments have been carried out for g 
somewhere in the range of panels (a) and (b). For (a), 
there is a very wide first integer SIRS on the upward 
sweep but none visible the downward direction. In (b) 
and (c), there are SIRS's in both directions, but wider 
on the upward sweep. In case (d), which we show for 
completeness but believe to correspond to an 



unattainable large coupling, there are no detectable steps 
but several discontinuities in the IV characteristic which 
are discussed below. The cavity energy E is calculated on 
the decreasing sweep. It shows a resonant enhancement 
even when the IV's (on this downward sweep) show no 
indication of a SIRS. [This enhancement is also visible on 
the upward sweep, which we have not shown.] In panel 
(a), shows a resonance at a current corresponding to a 
half-integer SIRS, but the IV characteristics themselves 
show no clear evidence of such a SIRS. In cases (b) and 
(c), we find that at these currents some fraction of the 
junctions have phase-locked onto the n — 1/2 step while 
the others are in the {Vj)r — state. Another notewor- 
thy feature is that as g increases, the integer steps in Fig. 
[To| (a) - (d) acquire a noticeable nonzero slope, and also 
become more and more rounded near their lower edge. 

In order to shed some light on the IV characteristics 
of Fig. |l^ (d), we have looked at the (Vi)r's across the 
individual junctions. Depending on I / Ic, all the {Vi)r's 
may be different, they may all be equal, or they fall into 
two or three groups. For certain I's, some of the {Vi)r's 
are nonzero while others vanish. This last behavior pre- 
sumably arises from the disorder in the critical currents. 



IV. DISCUSSION 

A. Comparison Between Calculated Results and 
Experiment 



w 



compare the present results to experi- 
Most of the published experiments thus far 



11 



Qj =100 



Qj=20 



q;=2 



Qj =0.10 



Qj =0.05 




0.7 1.0 
l/lo 
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(d) and (e), which correspond to overdamped junctions. The arrows indicate whether the trace is calculated for increasing or 
decreasing current. Lower panels show the time-averaged total energy E = (a|j + a|)r in the cavity, calculated as a function of 
decreasing current only. 



have been carried out on two-dimensional arrays. Their 
main features include the following: 

(a) . When the array is driven by a current, the IV 
characteristics show self-induced resonant steps. 

(b) . These steps are reported for any number of active 
junctions Na- 

(c) . Above a critical threshold number Nc of active 
junctions, the a. c. power output (i. e., the energy in 
the cavity) increases quadratically with Na- When the 
Na is increased through the threshold, the detected a. c. 
power in the cavity jumps by several orders of magnitude 
at the threshold. 

(d) . The array can be experimentally tuned so that 
different numbers of rows (i. e., different numbers of 
active junctions) are on the n = 1 SIRS. 

(e) . When Na junctions are on a SIRS and the current 
drive is varied, the Pac versus P^c curve is quadratic for 
low Pac and linear for high P^c- 

Our numerical results show all five of these features 
for a one dimensional array. Thus, they suggest that the 
behavior seen in the 2D experiments shoulcL-be visible 
even for a ID system. Indeed, a recent reporter suggests 
that all the features (a) - (e) are indeed experimentally 
observable in ID. 

We now elaborate on some of these points. The SIRS's 
emerge naturally from our equations of motion [Eqs. ( p^ ) 
and (pq)] . Another notable point is that we can numeri- 
cally control the number of active junctions Na by tuning 
the initial conditions. This tuning is possible because the 
junctions are underdamped and have an applied current 
regime within which they are bistable. The chosen Na 



determines whether the array is above or below the co- 
herence threshold Nc- If Na > Nc, then we usually find 
that, when the junctions lock onto a SIRS, they all lock 
onto the same, n = 1 step (first integer step). The voltage 
drop across the array is then {V)r / {N RIc) = Na^/Qj- 
Thus, the same array can produce an IV characteristic 
with multiple branches, each corresponding to a different 
number of SIRS's. This behavior is in agreement with 
the behavior seen in Ref. [17]. 

If Na < Nc, then our calculations still produce inte- 
ger SIRS's, but these steps are not coherent with one 
another. That is, although each junction is individu- 
ally locked onto the same fundamental frequency, which 
is close to the frequency fl of the cavity, the active junc- 
tions are out of phase with one another, and hence do not 
generate an energy in the cavity which varies quadrati- 
cally with Na- Also, even above the coherence threshold 
{Na > Nc), if the junctions are not locked on the steps, 
the array is not coherent at the coupling constant which 
produces the steps - that is, the power spectrum is remi- 
niscent of that of an array of independent junctions, and 
does not show a series of multiples of a single fundamen- 
tal frequency. Under these off-step conditions, the array 
can be made coherent, but only if the coupling constant 
is increased by several orders of magnitude above that 
needed to produce the SIRS's. 

Under some conditions, our calculations yield not only 
the first integer SIRS's but also overtone steps (higher 
integer steps), and fractional steps. The widths of our 
fractional steps are extremely small, and the steps are 
obtainable only by a delicate tuning of the current, initial 
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FIG. 10. This Figure illustrates the effects of changing the coupling parameter g while holding the other parameters fixed. 
Panels (a) - (d) correspond to g = 1.5 x 10"®, 4 x 10"^, 4 x 10"^, and 4 x 10"\ all with Qj = ^20, N = 10, and A = 0.05. 
Top panels: Time-averaged total voltage across the array, {V}t/{NRIc), versus current I/Ic- Arrows indicate the direction of 
current sweep. Bottom panels: total time-averaged energy E in the cavity as a function of I/Ic, all calculated for decreasing 
current bias. 



conditions, and current sweep rate. This sensitivity may 
explain why these fractional steps have not, as yet, been 
detected experimentally, though the overtone steps have 
been foundcj. 

Not only the general features but even some of the 
details of our calculations seem to agree well with exper- 
iment. For example, the results in Fig. ^(a) show the 
variation of a. c. power (that is, the electromagnetic en- 
ergy in the cavity) with the input d. c. power. The differ- 
ent curves correspond to distinct number of active junc- 
tions, Na, for this particular array. All the curves show 
a gradual, nearly parabolic onset but become nearly lin- 
ear at higher input power (that is, near the high-current 
edge of the step). The main difference between the cases 
Na > Nc and Na < Nc, is the behavior of the energy in 
the cavity after the SIRS becomes unstable (for increas- 
ing I/Q. When Na < N^, wc find that E - 10"^ at 
such input powers, while in the opposite case E 0.1. 
(This behavior is not shown in the Figure.) Very similar 
behavior to that shown in Fig. H(a) has recently been re- 
ported experimentally in Ref. [43], and is shown in Fig. 
^(b) . The similarity between the results of Rcf . [43] and 
the present work is apparent. A related experiment has 
also been reported in jsthich a 30% d. c. to a. c. conver- 
sion rate was achieveco. 



B. Qualitative Discussion of Underlying Physics 

We now briefly discuss the physics behind the present 
numerical results. First, the existence of a transition 



from incoherence to coherence, as a function of Na, re- 
sults from the "mean-field-like" nature of the interaction 
between the junctions and the cavity. Specifically, be- 
cause each junction is effectively coupled to every other 
junction via the cavity, the strength of the effective cou- 
pling increases with Na- Thus, for any g, a transition 
to coherence is to be expected for sufficiently large Na- 
A similar argument was made in the equilibrium case in 
Ref. [24]. 

Above the coherence transition, the self-induced reso- 
nant steps can also be qualitatively understood by refer- 
ring to the underlying equations ( ^ ) and (|4^). When 
a current is applied, it sets all the 7i's into motion, ac- 
cording to Eq. (^7|). If these 7i's all oscillate at the same 
fundamental frequency, they act as a driving term which 
causes an, and hence clr, to oscillate at the same fre- 
quency, according to Eq. (^). This then behaves like 
an a. c. current drive in Eq. (^). The combined d. c. 
and a. c. drives in Eq. (|4^) produce SIRS's, just as a com- 
bined d. c. and a. c. current produce Shapiro steps in a 
conventional Josephson junction. This same picture also 
makes it clear why the cavity energy increases quadrat- 
ically with Na above the threshold: in this regime, the 
"inhomogeneous" term on the right-hand side of Eq. ( |4^ ) 
is proportional to Na and, therefore, so is aj^. The whole 
process occurs self-consistently because the two equations 
are coupled. The effective "a. c. driving current" aji in 
Eq. (^) is also proportional to Na- Since the height of 
the first integer Shapiro step in a conventional junction 
is proportional to Ji{alac) where lac is the amplitude of 
the a. c. driving current and a is a constant related to the 
frequency, one might expect that the width of the SIRS's 
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would have an oscillatory dependence on Na- There are 
some slight hints of this behavior in our numerical results 
[cf. Fig. 1(a)]. 

This description also suggests why the steps occur even 
in one-dimensional arrays. Their occurrence depends, 
not on the dimensionality of the array, but only on the 
existence of a suitable induced a. c. drive. Indeed, such 
steps have recently been reported in ID arraysLJ, consis- 
tent with the present model. The in-plane magnetic field 
used in the earlier experiments is apparently needed only 
to lower the Josephson critical currents sufficiently that 
the resonant frequency SI occurs in the bistable region of 
the IV characteristics. 

All the numerical results in the present paper are ob- 
tained in the "semi-classical" regime, where the various 
operators are regarded as c- numbers. It would be of inter- 
est to study the array dynamics of the array in the quan- 
tum regime, where the number of photons is small. A 
recent numerical study of this kind (but only for the equi- 
librium properties) has been carried out for a SQUID iip_a 
resonant cavity (without resistively-shunted damping)Ej. 

In summary, we have derived the Heisenberg equations 
of motion for a model Hamiltonian which describes a one- 
dimensional array of underdamped Josephson junctions 
coupled to a resonant cavity. We have numerically solved 
these equations in the classical limit, valid in the limit of 
large numbers of photons in the cavity. In the presence 
of a d. c. current drive, we find numerically that (i) the 
array exhibits self- induced resonant steps (SIRS), simi- 
lar to Shapiro steps in conventional arrays; (ii) there is 
a transition between an unsynchronized and a synchro- 
nized state as the number of active junctions is increased 
while other parameters are held fixed; and (iii) when the 
array is biased on the first integer SIRS, the total energy 
increases quadratically with number of active junctions. 
Our results are in quite detailed agreement with exper- 
iment, even though the experiments are largely carried 
out in 2D. Thus, the present model strongly suggests 
that a 2D array is not necessary in order to obtain the 
observed SIRS's. The results also strongly suggest that 
the experimental data considered here can be understood 
in terms of a model involving strictly classical equations 
of motion, without the necessity of introducing new, non- 
classical physics. 
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